
## PEASANT RESISTANCE AND ECONOMIC AFFLUENCE: LESSONS FROM PARAGUAY
## Appendix C: Evidence on the Mechanisms
## January 2024

## Code Authors:
## Germán Feierherd (gfeierherd@udesa.edu.ar) 
## Jorge Mangonnet (jorge.mangonnet@vanderbilt.edu)


#### Preamble ------------------------------------------------------------------

## Workspace setup
rm(list = ls())		         # clear all objects in memory
options(max.print = 5000)  # expand display
options(scipen = 999)      # disable sci notation
dev.off()                  # reload graphic device
cat("\014")                # clear console

## Load/install packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl,
  interflex,
  marginaleffects,
  haven,
  readr,
  reshape2,
  magrittr,
  dplyr,
  plm,
  lfe,
  pcse,
  lmtest,
  stargazer,
  xtable,
  texreg,
  estimatr, 
  fixest,
  margins, 
  sjPlot, 
  multiwayvcov,
  #  rgdal,
  ggplot2,
  wesanderson, 
  fixest,
  did,
  sf)

## Set working directory
gf <- "~/Dropbox/Working_Papers/investigacion protesta Paraguay/"; main <- gf  # german's path
jm <- "~/Dropbox/Research/investigacion protesta Paraguay/"; main <- jm        # jorge's path
setwd(main)

## Functions
source(paste0(main, "4_code/functions_paraguay.R"))

## Load full dataset
load(paste0(main, "4_code/DataSetPY.Rda"))


#### Appendix C: Evidence on the Mechanisms ------------------------------------

# Municipality-level data frame
dta_mun <- dta %>% dplyr::select(codigo,departamento,
                                 ln_suitability, ln_soy_suit, 
                                 ln_mze_suit, ln_sug_suit,
                                 ln_distance_asuncion,ln_distance_hubs,
                                 ln_total_conflicts,
                                 distance_asuncion,distance_hubs,
                                 matches("planted|produced")) %>% distinct()
dta_mun[is.na(dta_mun)] <- 0


#### Table C.1: Planted Hectares and Land Suitability ----

# Model 1: Total Planted (2008-1991)
m1 <- lm(total_planted_diff ~ ln_suitability + factor(departamento), data = dta_mun); summary(m1)
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC0")); se1

# Model 2: Soybean Planted (2008-1991)
m2 <- lm(soybean_planted_diff ~ ln_soy_suit + factor(departamento), data = dta_mun); summary(m2)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC0")); se2

# Model 3: Maize Planted (2008-1991)
m3 <- lm(maize_planted_diff ~ ln_mze_suit + factor(departamento), data = dta_mun); summary(m3)
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC0")); se2

# Model 4: Sugar Planted (2008-1991)
m4 <- lm(sugar_planted_diff ~ ln_sug_suit + factor(departamento), data = dta_mun); summary(m4)
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC0")); se4

# Model 5: Soybean Planted (2014-1991)
m5 <- lm(soybean_planted_diff_14 ~ ln_soy_suit + factor(departamento), data = dta_mun); summary(m5)
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC0")); se5

# Create table
stargazer(
  m1, m2, m3, m4, m5,
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4]),  
  covariate.labels = c("Suitability",
                       "Soybean Suitability",
                       "Maize Suitability",
                       "Sugar Suitability"), 
  add.lines = list(c("","","","","")), 
  #c("Department FE", rep("Yes", 5)),
  #c("Year FE", rep("No", 5))),
  omit = c("factor"),
  omit.stat = c("adj.rsq", "ser", "f", "lr", "aic"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  column.labels = c("\\text{\\normalfont Total (08-91)}",
                    "\\text{\\normalfont Soybeans (08-91)}",
                    "\\text{\\normalfont Maize (08-91)}",
                    "\\text{\\normalfont Sugar (08-91)}",
                    "\\text{\\normalfont Soybeans (14-91)}"),
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  font.size = "scriptsize",
  #table.placement = "p!",
  style = "ajps", 
  title = "Percentage Change in Planted Hectares and Land Suitability", 
  label = "planted",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{13.3cm}{\\textit{Note}: All models include departmental effects. Huber-White robust standard errors in parentheses. The data come from the 2008 and 1991 agricultural censuses and CAPECO.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_planted_dv2.tex")
)


#### Table C.2: Peasant Resistance and Planted Hectares ----

# Model 1: Total Planted (2008-1991)
m1 <- lm(ln_total_conflicts ~ total_planted_diff + factor(departamento), data = dta_mun); summary(m1)
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC0")); se1

# Model 2: Soybean Planted (2008-1991)
m2 <- lm(ln_total_conflicts ~ soybean_planted_diff + factor(departamento), data = dta_mun); summary(m2)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC0")); se2

# Model 3: Maize Planted (2008-1991)
m3 <- lm(ln_total_conflicts ~ maize_planted_diff + factor(departamento), data = dta_mun); summary(m3)
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC0")); se2

# Model 4: Sugar Planted (2008-1991)
m4 <- lm(ln_total_conflicts ~ sugar_planted_diff + factor(departamento), data = dta_mun); summary(m4)
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC0")); se4

# Model 5: Soybean Planted (2014-1991)
m5 <- lm(ln_total_conflicts ~ soybean_planted_diff_14 + factor(departamento), data = dta_mun); summary(m5)
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC0")); se5

# Create table
stargazer(
  m1, m2, m3, m4, m5,
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4]),  
  covariate.labels = c("Total (08-91)",
                       "Soybean (08-91)",
                       "Maize (08-91)",
                       "Sugar (08-91)",
                       "Soybean (14-91"), 
  add.lines = list(c("","","","","")), 
  #c("Department FE", rep("Yes", 5)),
  #c("Year FE", rep("No", 5))),
  omit = c("factor"),
  omit.stat = c("adj.rsq", "ser", "f", "lr", "aic"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  font.size = "small",
  #table.placement = "p!",
  style = "ajps", 
  title = "Sum of Peasant Resistance and Percentage Change in Planted Hectares", 
  label = "resist",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{12.5cm}{\\textit{Note}: All models include departmental effects. Huber-White robust standard errors in parentheses. The data come from the 2008 and 1991 agricultural censuses and CAPECO.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_planted_iv2.tex")
)


#### Figure C.1: Local Regressions of Peasant Resistance, Planted Hectares, and Land Suitability ----

# Plot (a): % planted hectares (2008-1991)
lwr.planted <- ggplot(dta_mun, aes(ln_suitability, total_planted_diff)) + 
  geom_smooth(method = "lm", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Land Suitability") + ylab(expression(Delta * "% Planted Hectares (2008-1991)"))

pdf(file = paste0(main, "1_writing/2023/figures/loess_planted.pdf"))
lwr.planted
dev.off()

# Plot (b): sum of peasant resistance
lwr.resist <- ggplot(dta_mun, aes(total_planted_diff, ln_total_conflicts)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab(expression(Delta * "% Planted Hectares (2008-1991)")) + ylab("Sum of Resistance Events (IHS)")

pdf(file = paste0(main, "1_writing/2023/figures/loess_resist.pdf"))
lwr.resist
dev.off()


#### Figure C.2: Local Regressions of Planted Hectares and Land Suitability by Capital-Intensive Crop ----

# Plot (a): soybean (2008-1991)
lwr.soy <- ggplot(dta_mun, aes(ln_soy_suit, soybean_planted_diff)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Soybean Suitability") + ylab(expression(Delta * "% 2008-1991"))

pdf(file = paste0(main, "1_writing/2023/figures/loess_soy_planted.pdf"))
lwr.soy
dev.off()

# Plot (b): soybean (2014-1991)
lwr.soy14 <- ggplot(dta_mun, aes(ln_soy_suit, soybean_planted_diff_14)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Soybean Suitability") + ylab(expression(Delta * "% 2014-1991"))

pdf(file = paste0(main, "1_writing/2023/figures/loess_soy14_planted.pdf"))
lwr.soy14
dev.off()

# Plot (c): maize (2008-1991)
lwr.mze <- ggplot(dta_mun, aes(ln_mze_suit, maize_planted_diff)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Maize Suitability") + ylab(expression(Delta * "% 2008-1991"))

pdf(file = paste0(main, "1_writing/2023/figures/loess_mze_planted.pdf"))
lwr.mze
dev.off()

# Plot (d): sugar (2008-1991)
lwr.sug <- ggplot(dta_mun, aes(ln_sug_suit, sugar_planted_diff)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  coord_cartesian(xlim = c(-4.2,-4.096)) +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Sugar Suitability") + ylab(expression(Delta * "% 2008-1991"))

pdf(file = paste0(main, "1_writing/2023/figures/loess_sug_planted.pdf"))
lwr.sug
dev.off()


#### Figure C.3: Local Regressions of Peasant Resistance and Planted Hectares by Capital-Intensive Crop ----

# Plot (a): soybean (2008-1991)
lwr.soy <- ggplot(dta_mun, aes(soybean_planted_diff, ln_total_conflicts)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab(expression(Delta * "% Soybean Hectares 2008-1991")) + ylab("Sum of Resistance Events")

pdf(file = paste0(main, "1_writing/2023/figures/loess_soy_resist.pdf"))
lwr.soy
dev.off()

# Plot (b): soybean (2014-1991)
lwr.soy14 <- ggplot(dta_mun, aes(soybean_planted_diff_14, ln_total_conflicts)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab(expression(Delta * "% Soybean Hectares 2014-1991")) + ylab("Sum of Resistance Events")

pdf(file = paste0(main, "1_writing/2023/figures/loess_soy14_resist.pdf"))
lwr.soy14
dev.off()

# Plot (c): maize (2008-1991)
lwr.mze <- ggplot(dta_mun, aes(maize_planted_diff, ln_total_conflicts)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab(expression(Delta * "% Maize Hectares 2008-1991")) + ylab("Sum of Resistance Events")

pdf(file = paste0(main, "1_writing/2023/figures/loess_mze_resist.pdf"))
lwr.mze
dev.off()

# Plot (d): sugar (2008-1991)
lwr.sug <- ggplot(dta_mun, aes(sugar_planted_diff, ln_total_conflicts)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  #coord_cartesian(xlim = c(-4.2,-4.096)) +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab(expression(Delta * "% Sugar Hectares 2008-1991")) + ylab("Sum of Resistance Events")

pdf(file = paste0(main, "1_writing/2023/figures/loess_sug_resist.pdf"))
lwr.sug
dev.off()


#### Table C.3: Peasant Resistance to Income Reductions ----

# Model 1: Prices x Suitability, w/o controls
m1 <- plm(ln_conflict_inc ~ I(ln_suitability*ln_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability, w/ controls
m2 <- plm(ln_conflict_inc ~ I(ln_suitability*ln_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_conflict_inc ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_conflict_inc ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_conflict_inc ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_conflict_inc ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices $\\times$ Suitability", 
                       "Prices $\\times$ Settlements", 
                       "Prices $\\times$ Suitability $\\times$ Settlements", 
                       "Prices $\\times$ Committees", 
                       "Prices $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
  #c("Municipality FE", rep("Yes", 6)),
  #c("Year FE", rep("Yes", 6))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps", 
  title = "Peasant Resistance to Income Reductions", 
  label = "income",
  notes.append = FALSE,
  notes = c("\\parbox[t]{12.5cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. The dependent variable are all the events of peasant resistance against reductions in rural income (e.g, low wages or insufficient rural credit). The data come from the authors' archival database.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_income.tex")
)


#### Figure C.4: Marginal Effect of Agricultural Prices on Peasant Resistance to Income Reductions by Land Suitability ----

# Based on model 1
m1 <- plm(ln_conflict_inc ~ ln_suitability*ln_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_income.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price", moderator="ln_suitability", interaction="ln_suitability:ln_price", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure C.5: Marginal Effect of Agricultural Prices on Peasant Resistance to Income Reductions by Land Suitability, Subsistence Agriculture, and Organizational Capacities ----

# Based on model 3 & 5
m3 <- plm(ln_conflict_inc ~ ln_suitability*ln_price + ln_price*subsistence_farms_median + ln_suitability*ln_price*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_conflict_inc ~ ln_suitability*ln_price + ln_price*orgs_locales_median + ln_suitability*ln_price*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_income.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_income.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agriculutral Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6, covMat1, covMat3, covMat5)


#### Table C.4: Peasant Resistance to Indigenous Peoples ----

# Model 1: Prices x Suitability, w/o controls
m1 <- plm(ln_peasvsindig ~ I(ln_suitability*ln_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability, w/ controls
m2 <- plm(ln_peasvsindig ~ I(ln_suitability*ln_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_peasvsindig ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_peasvsindig ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_peasvsindig ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_peasvsindig ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices $\\times$ Suitability", 
                       "Prices $\\times$ Settlements", 
                       "Prices $\\times$ Suitability $\\times$ Settlements", 
                       "Prices $\\times$ Committees", 
                       "Prices $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
  #c("Municipality FE", rep("Yes", 6)),
  #c("Year FE", rep("Yes", 6))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps", 
  title = "Peasant Resistance to Indigenous Peoples", 
  label = "peasvsindig",
  notes.append = FALSE,
  notes = c("\\parbox[t]{14cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. The dependent variable are all the events of peasant resistance to reductions in rural income (e.g, low wages or insufficient rural credit). The data come from the authors' archival database.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_peasvsindig.tex")
)


#### Figure C.6: Marginal Effect of Agricultural Prices on Peasant Resistance to Indigenous Peoples by Land Suitability ----

# Based on model 1
m1 <- plm(ln_peasvsindig ~ ln_suitability*ln_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_peasvsindig.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price", moderator="ln_suitability", interaction="ln_suitability:ln_price", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure C.6: Marginal Effect of Agricultural Prices on Peasant Resistance to Indigenous Peoples by Land Suitability, Subsistence Agriculture, and Organizational Capacities ----

# Based on model 3 & 5
m3 <- plm(ln_peasvsindig ~ ln_suitability*ln_price + ln_price*subsistence_farms_median + ln_suitability*ln_price*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_peasvsindig ~ ln_suitability*ln_price + ln_price*orgs_locales_median + ln_suitability*ln_price*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_peasvsindig.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_peasvsindig.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agriculutral Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6, covMat1, covMat3, covMat5)


#### Table C.5: Landowner Protest Against the Government ----

# Model 1: Prices x Suitability, w/o controls
m1 <- plm(ln_landvsstate ~ I(ln_suitability*ln_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability, w/ controls
m2 <- plm(ln_landvsstate ~ I(ln_suitability*ln_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_landvsstate ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_landvsstate ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4)
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_landvsstate ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_landvsstate ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices $\\times$ Suitability", 
                       "Prices $\\times$ Settlements", 
                       "Prices $\\times$ Suitability $\\times$ Settlements", 
                       "Prices $\\times$ Committees", 
                       "Prices $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
  #c("Municipality FE", rep("Yes", 6)),
  #c("Year FE", rep("Yes", 6))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps", 
  title = "Landowner Protest Against the State", 
  label = "landvsstate",
  notes.append = FALSE,
  notes = c("\\parbox[t]{12.75cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. The dependent variable are all the events of landowner protest against and the government (e.g., tractor roadblocks demanding lower taxes or lower gas prices). The data come from the authors' archival database.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_landvsstate.tex")
)


#### Figure C.8: Marginal Effect of Agricultural Prices on Landowner Protest Against the Government by Land Suitability ----

# Based on model 1
m1 <- plm(ln_landvsstate ~ ln_suitability*ln_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_landvsstate.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price", moderator="ln_suitability", interaction="ln_suitability:ln_price", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure C.9: Marginal Effect of Agricultural Prices on Landowner Protest Against the Government by Land Suitability, Subsistence Agriculture, and Organizational Capacities ----

# Based on model 3 & 5
m3 <- plm(ln_landvsstate ~ ln_suitability*ln_price + ln_price*subsistence_farms_median + ln_suitability*ln_price*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_landvsstate ~ ln_suitability*ln_price + ln_price*orgs_locales_median + ln_suitability*ln_price*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_landvsstate.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_landvsstate.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agriculutral Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6, covMat1, covMat3, covMat5)






#### Table C.6: Peasant Resistance and Labor-intensive Agriculture ----

# Model 1: Prices x Suitability, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability_pbo*ln_price_pbo) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability_pbo*ln_price_pbo) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_conflict ~ I(ln_suitability_pbo*ln_price_pbo) + I(ln_price_pbo*subsistence_farms_median) + I(ln_suitability_pbo*ln_price_pbo*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_conflict ~ I(ln_suitability_pbo*ln_price_pbo) + I(ln_price_pbo*subsistence_farms_median) + I(ln_suitability_pbo*ln_price_pbo*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_conflict ~ I(ln_suitability_pbo*ln_price_pbo) + I(ln_price_pbo*orgs_locales_median) + I(ln_suitability_pbo*ln_price_pbo*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_conflict ~ I(ln_suitability_pbo*ln_price_pbo) + I(ln_price_pbo*orgs_locales_median) + I(ln_suitability_pbo*ln_price_pbo*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices $\\times$ Suitability", 
                       "Prices $\\times$ Settlements", 
                       "Prices $\\times$ Suitability $\\times$ Settlements", 
                       "Prices $\\times$ Committees", 
                       "Prices $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
  #c("Municipality FE", rep("Yes", 6)),
  #c("Year FE", rep("Yes", 6))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps",  
  title = "Peasant Resistance and Labor-Intensive Agriculture", 
  label = "noflex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{13.25cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Prices is the log of Paraguay's average producer price of cotton, tobacco, and yerba mate (Paraguay's labor-intensive crops) in U.S. dollars per metric ton. Suitability is the log of Paraguay's average potential yield for cotton, tobacco, and yerba mate in metric tons per hectare. The data come from FAOSTAT and GAEZ, respectively.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_noflex.tex")
)


#### Figure C.10: Marginal Effect of Agricultural Prices on Peasant Resistance to Labor-Intensive Agriculture by Land Suitability ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability_pbo*ln_price_pbo + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_noflex.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price_pbo", moderator="ln_suitability_pbo", interaction="ln_suitability_pbo:ln_price_pbo", varcov=covMat1, minimum=-1.7, maximum=-1, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure C.11: Marginal Effect of Agricultural Prices on Peasant Resistance to Labor-Intensive Agriculture by Land Suitability, Subsistence Agriculture, and Organizational Capacities ----

# Based on model 3 & 5
m3 <- plm(ln_conflict ~ ln_suitability_pbo*ln_price_pbo + ln_price_pbo*subsistence_farms_median + ln_suitability_pbo*ln_price_pbo*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_conflict ~ ln_suitability_pbo*ln_price_pbo + ln_price_pbo*orgs_locales_median + ln_suitability_pbo*ln_price_pbo*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_noflex.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price_pbo", moderator1="ln_suitability_pbo", moderator2="subsistence_farms_median", varcov=covMat3, minimum=-1.7, maximum=-1, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_noflex.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price_pbo", moderator1="ln_suitability_pbo", moderator2="orgs_locales_median", varcov=covMat5, minimum=-1.7, maximum=-1, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agriculutral Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6, covMat1, covMat3, covMat5)


#### Table C.7: Peasant Resistance and Food Inflation ----

# Model 1: Inflation x Suitability, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*inflation) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Inflation x Suitability, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*inflation) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Inflation x Suitability x Settlements, w/o controls
m3 <- plm(ln_conflict ~ I(ln_suitability*inflation) + I(inflation*subsistence_farms_median) + I(ln_suitability*inflation*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Inflation x Suitability x Settlements, w/ controls
m4 <- plm(ln_conflict ~ I(ln_suitability*inflation) + I(inflation*subsistence_farms_median) + I(ln_suitability*inflation*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Inflation x Suitability x Committees, w/o controls
m5 <- plm(ln_conflict ~ I(ln_suitability*inflation) + I(inflation*orgs_locales_median) + I(ln_suitability*inflation*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", na.action = na.omit, data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Inflation x Suitability x Committees, w/ controls
m6 <- plm(ln_conflict ~ I(ln_suitability*inflation) + I(inflation*orgs_locales_median) + I(ln_suitability*inflation*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", na.action = na.omit, data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Consumer Price Index $\\times$ Suitability", 
                       "Consumer Price Index $\\times$ Settlements", 
                       "Consumer Price Index $\\times$ Suitability $\\times$ Settlements", 
                       "Consumer Price Index $\\times$ Committees", 
                       "Consumer Price Index $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")), 
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps", 
  title = "Peasant Resistance and Food Inflation", 
  label = "inflation",
  notes.append = FALSE,
  notes = c("\\parbox[t]{14cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Food inflation is the annual change in the cost of a basic food basket. The data for food inflation come from the World Bank.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_inflation.tex")
)


#### Figure C.12: Marginal Effect of Food Inflation on Peasant Resistance by Land Suitability ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability*inflation + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_inflation.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="inflation", moderator="ln_suitability", interaction="ln_suitability:inflation", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Food Inflation")
dev.off()


#### Figure C.13: Marginal Effect of Food Inflation on Peasant Resistance by Land Suitability, Subsistence Settlements, and Peasant Committees ----

# Based on model 3 & 5
m3 <- plm(ln_conflict ~ ln_suitability*inflation + inflation*subsistence_farms_median + ln_suitability*inflation*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_conflict ~ ln_suitability*inflation + inflation*orgs_locales_median + ln_suitability*inflation*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_inflat.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="inflation", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Food Inflation")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_inflat.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="inflation", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Food Inflation")
dev.off()


#### Table C.8: Peasant Commercial Agriculture ----

# Model 1: Prices x Suitability x Cotton Farms, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*cotton_farms20_median) + I(ln_suitability*ln_price*cotton_farms20_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability x Cotton Farms, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*cotton_farms20_median) + I(ln_suitability*ln_price*cotton_farms20_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = "codigo", model = "within", effect = "individual", data = dta, na.action = na.exclude); summary(m2)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Create table 
stargazer(
  m1, m2, 
  se = list(se1[,2], se2[,2]), 
  p = list(se1[,4], se2[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices $\\times$ Land Suitability", 
                       "Prices $\\times$ Peasant Cotton Farms", 
                       "Prices $\\times$ Land Suitability $\\times$ Peasant Cotton Farms"),
  add.lines = list(c("",""), 
                   c("Controls", "No", "Yes")),
  #c("Municipality FE", rep("Yes", 2)),
  #c("Year FE", rep("Yes", 2))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  float = TRUE, 
  #float.env = "sidewaystable",
  style = "ajps", 
  title = "Peasant Commercial Agriculture", 
  label = "cotton",
  notes.append = FALSE,
  notes = c("\\parbox[t]{16cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Peasant cotton farms is a dummy variable measuring whether the number of small farms (5 hectares or less) planting cotton is greater than the median. The data for cotton farms come from the 1991 agricultural census.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_cotton.tex")
)


#### Figure C.14: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability and Peasant Commercial Agriculture ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*cotton_farms20_median + ln_suitability*ln_price*cotton_farms20_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 

# Covariance matrices for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_cotton.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m1, effect="ln_price", moderator1="ln_suitability", moderator2="cotton_farms20_median", varcov=covMat1, minimum=.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()

remove(m1, m2, se1, se2, covMat1)


#### Table C.9: IBR Colonization ----

#### IBR peasant colonies ----

# Model 1: Prices x Suitability x IBR Colonies, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*old_colonies_median) + I(ln_suitability*ln_price*old_colonies_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability x IBR Colonies, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*old_colonies_median) + I(ln_suitability*ln_price*old_colonies_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = "codigo", model = "within", effect = "individual", data = dta); summary(m2)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Create table
stargazer(
  m1, m2,
  se = list(se1[,2], se2[,2]), 
  p = list(se1[,4], se2[,4]),  
  covariate.labels = c("Prices $\\times$ Land Suitability",
                       "Prices $\\times$ IBR Colonies",
                       "Prices $\\times$ Land Suitability $\\times$ IBR Colonies"),
  add.lines = list(c("",""), 
                   c("Controls", c("No", "Yes"))),
  #c("Municipality FE", rep("Yes", 2)),
  #c("Year FE", "Yes", rep("Yes", 2))),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  #table.placement = "p!",
  style = "ajps", 
  title = "IBR Colonization", 
  label = "colonies",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{14cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. IBR colonies is a dummy variable measuring whether the number of peasant colonies founded by the IBR in 1963-1989 is greater than the median. The data come from \\citet{RojasAreco2017}.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_colonies.tex")
)


#### Figure C.15: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability and IBR Colonization ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*old_colonies_median + ln_suitability*ln_price*old_colonies_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_colonies.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m1, effect="ln_price", moderator1="ln_suitability", moderator2="old_colonies_median", varcov=covMat1, minimum=.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()

remove(m1, m2, se1, se2, covMat1)


#### Table C.10: Landowner Associations ----

# Model 1: Prices x Suitability x Landowner Associations, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*landowner_orgs_median) + I(ln_suitability*ln_price*landowner_orgs_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability x Landowner Associations, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*landowner_orgs_median) + I(ln_suitability*ln_price*landowner_orgs_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = "codigo", model = "within", effect = "individual", data = dta); summary(m2)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Create table
stargazer(
  m1, m2,
  se = list(se1[,2], se2[,2]), 
  p = list(se1[,4], se2[,4]),  
  covariate.labels = c("Prices $\\times$ Land Suitability",
                       "Prices $\\times$ Landowner Association",
                       "Prices $\\times$ Land Suitability $\\times$ Landowner Association"),
  add.lines = list(c("",""), 
                   c("Controls", c("No", "Yes"))),
  #c("Municipality FE", rep("Yes", 2)),
  #c("Year FE", "Yes", rep("Yes", 2))),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  #table.placement = "p!",
  style = "ajps", 
  title = "Landowner Associations", 
  label = "landowner",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{16cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Landowner association is a dummy variable measuring whether the number of landowner associations affiliated with the ARP or UGP is greater than the median. The data come from the ARP and UGP's list of affiliates and district offices.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_landowner.tex")
)


#### Figure C.16: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability and Landowner Associations ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*landowner_orgs_median + ln_suitability*ln_price*landowner_orgs_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_landowner.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m1, effect="ln_price", moderator1="ln_suitability", moderator2="landowner_orgs_median", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()

remove(m1, m2, se1, se2, covMat1)




